% This code: transform shapefile with road network into a graph

function [ road_network,unique_nodes,unique_nodes_nat ] = transform_network_to_graph( code_control,roads,country )

    
    %% return all the coordinates defined in the road network shapefile
    all_coord  = [];  % coordinates of each node in road network
    for j=1:length( roads )

        [ j length( roads ) ]
        
        temp = [ roads(j).X( 1:length( roads(j).X )-1 )' ...
                 roads(j).Y( 1:length( roads(j).X )-1 )' ];
        all_coord = [ all_coord;temp ];            
    end

    % identify the unique nodes
    [ unique_nodes,~,ic ] = unique( all_coord,'rows' );  
           % unique_nodes = coordinates of the unique nodes. nodes are numbered by its position in this vector
           % all_coord = unique_nodes( ic ), so ic indicates the unique node # associated with each elemen of all_coord
           % unique_nodes = all_coord( ia )   

    % assign an integer node index to each coordinate
    k=1;
    for j=1:length(roads)

        roads(j).node_index = ic( k:k+length( roads(j).X )-2 );   % these are the nodes indexes of all the nodes in roads(j)
        k = k+length( roads(j).X )-1;
        
    end
    
    %% return all the coordinates corresponding to national nodes defined in the road network shapefile
    
    all_coord_nat = [];  % coordinates of each node in road network    
    roads_nat = roads( logical( cell2mat( {roads.use} ) ) );   % national roads only
        
    for j=1:length( roads_nat )
        
        [ j length( roads_nat ) ]

        temp = [ roads_nat(j).X( 1:length( roads_nat(j).X )-1 )' ...
                 roads_nat(j).Y( 1:length( roads_nat(j).X )-1 )' ];              
             
        all_coord_nat = [ all_coord_nat;temp ];   
    end

    % identify the unique nodes
    [ unique_nodes_nat,~,ic_nat ] = unique( all_coord_nat,'rows' );  
           % unique_nodes = coordinates of the unique nodes. nodes are numbered by its position in this vector
           % all_coord = unique_nodes( ic ), so ic indicates the unique node # associated with each elemen of all_coord
           % unique_nodes = all_coord( ia )   

    % assign an integer node index to each coordinate
    k=1;
    for j=1:length( roads_nat )
        
        [ j length( roads_nat ) ]

        roads_nat(j).node_index = ic_nat( k:k+length( roads_nat(j).X )-2 );   % these are the nodes indexes of all the nodes in roads(j)
        k = k+length( roads_nat(j).X )-1;
        
    end

    
     
    %%
    %{
    note: 

    [ unique_nodes(roads(j).node_index,:); nan nan ] = [roads(j).X' roads(j).Y']

    

    unique_nodes( roads(j).node_index,: ) = [ roads(j).X(1:length( roads(j).X )-1 )' ...
                                              roads(j).Y(1:length( roads(j).X )-1 )']
    %}

    % define graph
    road_network_s = [];
    road_network_t = [];
    edge_lanes = [];
    edge_use = [];
    edge_distance = [];
    edge_weights = [];
    
    for j=1:length(roads)
        
        [ j length(roads) ]

        road_network_s = [ road_network_s;roads(j).node_index( 1:length( roads(j).node_index )-1 ) ];   % starting node in each edge
        
        road_network_t = [ road_network_t;roads(j).node_index( 2:length( roads(j).node_index ) ) ];     % final node in each edge
       
        edge_lanes = [ edge_lanes;ones( length( roads(j).node_index )-1,1 )*roads(j).lanes  ];  % number of actual lanes in each edge
        
        edge_use = [ edge_use;ones( length( roads(j).node_index )-1,1 )*roads(j).use  ];        % length (in KM) of each edge
        
        edge_distance = [ edge_distance;roads(j).distances' ];
        
        edge_weights = [  edge_weights;roads(j).weights' ];
            
    end


%% graph with road network

EdgeTable = table( [road_network_s road_network_t],...   % nodes
                    edge_weights,...                     % weights
                    edge_distance,...                    % actual distance
                    edge_lanes,...              % lanes
                    edge_use,...    % use
                    'VariableNames',{'EndNodes','Weight','Distance','lanes','use'} );   % variable names
                
road_network = graph( EdgeTable );

% Fill in missing values with the mean weight, and missing lanes with one lane     
road_network.Edges.lanes(ismissing(road_network.Edges.lanes)) = ones( sum(ismissing(road_network.Edges.lanes)), 1);
road_network.Edges.Weight(ismissing(road_network.Edges.Weight)) = mean(road_network.Edges.Weight(~ismissing(road_network.Edges.Weight)));




